%----------------------------------------------
% Formatting and file path
%----------------------------------------------

clear;
clc;
format bank;

for j=70:100

    %----------------------------------------------
    % Simulate GE price effects by reducing A
    %----------------------------------------------

    load('../../Data/structural/output/Structural_Results_at_ss_assets_minus_land_of_rich.mat')
    A=A.*(1-(j/100));

    %----------------------------------------------
    % Estimate h*, l*, h' at rich ss level of ALL ASSETS MINUS LAND
    %----------------------------------------------

    k=43701;
    fk=a*(k^2) + b*k;

    psi_l_sim = ratio_psi_l.*psi_l;
    psi_h_sim = ratio_psi_l.*psi_h;

    l_opt = zeros(N,1);
    h_opt = zeros(N,1);
    hiredin_opt = zeros(N,1);
    profit_opt = zeros(N,1);

    for i=1:N
        profit = @(x) -A(i)*fk*(x(1) + x(3))^beta - w_bl(i)*x(2) + phi*w_bl(i)*x(3) + 0.5*(sqrt(psi_l_sim(i))*x(1) + sqrt(psi_h_sim(i))*x(2))^2;
        [argmax, max] = fmincon(profit, [0,0,0], [1,1,0], R, [], [], [0,0,0], [R,H(i),Hire_in_max]);
        l_opt(i) = argmax(1);
        h_opt(i) = argmax(2);
        hiredin_opt(i) = argmax(3);
        profit_opt(i) = -max;
    end

    clear max argmax

    % Check if it would be optimal to liquidate capital:

    for i=1:N
        profitalt = @(x) -rho*k -w_bl(i)*x + 0.5*psi_h_sim(i)*x^2;
        [argmax, max] = fmincon(profitalt, 0, [], [], [], [], 0, H(i));
        if -max > profit_opt(i)
            l_opt(i) = 0;
            hiredin_opt(i) = 0;
            h_opt(i) = argmax;
            profit_opt(i) = -max;
        end
    end

    clear max argmax

    % Trim :
    l_opt(l_opt < 1) = 0;
    h_opt(h_opt < 1) = 0;
    hiredin_opt(hiredin_opt < 1) = 0;

    % Optimal cases:
    case_opt = zeros(N,1);

    for i=1:N
        if l_opt(i) > 0 && h_opt(i) > 0 && hiredin_opt(i) > 0
            case_opt(i) = 1;
        elseif l_opt(i) > 0 && h_opt(i) > 0 && hiredin_opt(i) == 0
            case_opt(i) = 2;
        elseif l_opt(i) > 0 && h_opt(i) == 0 && hiredin_opt(i) > 0
            case_opt(i) = 3;
        elseif l_opt(i) > 0 && h_opt(i) == 0 && hiredin_opt(i) == 0
            case_opt(i) = 4;
        elseif l_opt(i) == 0 && h_opt(i) > 0 && hiredin_opt(i) > 0
            case_opt(i) = 5;
        elseif l_opt(i) == 0 && h_opt(i) > 0 && hiredin_opt(i) == 0
            case_opt(i) = 6;
        elseif l_opt(i) == 0 && h_opt(i) == 0
            case_opt(i) = 7;
        end
    end

    %----------------------------------------------
    % Estimate h*, l*, h' at baseline capital level
    %----------------------------------------------

    l_opt_bl = zeros(N,1);
    h_opt_bl = zeros(N,1);
    hiredin_opt_bl = zeros(N,1);
    profit_opt_bl = zeros(N,1);

    for i=1:N
        profit = @(x) -A(i)*fk_bl(i)*(x(1) + x(3))^beta - w_bl(i)*x(2) + phi*w_bl(i)*x(3) + 0.5*(sqrt(psi_l(i))*x(1) + sqrt(psi_h(i))*x(2))^2;
        [argmax, max] = fmincon(profit, [0,0,0], [1,1,0], R, [], [], [0,0,0], [R,H(i),Hire_in_max]);
        l_opt_bl(i) = argmax(1);
        h_opt_bl(i) = argmax(2);
        hiredin_opt_bl(i) = argmax(3);
        profit_opt_bl(i) = -max;
    end

    clear max argmax

    % Check if it would be optimal to liquidate capital:

    for i=1:N
        profitalt = @(x) -rho*k_bl_pre_trans(i) -w_bl(i)*x + 0.5*psi_h(i)*x^2;
        [argmax, max] = fmincon(profitalt, 0, [], [], [], [], 0, H(i));
        if -max > profit_opt_bl(i)
            l_opt_bl(i) = 0;
            hiredin_opt_bl(i) = 0;
            h_opt_bl(i) = argmax;
            profit_opt_bl(i) = -max;
        end
    end

    clear max argmax

    % Trim :
    l_opt_bl(l_opt_bl < 1) = 0;
    h_opt_bl(h_opt_bl < 1) = 0;
    hiredin_opt_bl(hiredin_opt_bl < 1) = 0;

    % Optimal cases:
    case_opt_bl = zeros(N,1);

    for i=1:N
        if l_opt_bl(i) > 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) > 0
            case_opt_bl(i) = 1;
        elseif l_opt_bl(i) > 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) == 0
            case_opt_bl(i) = 2;
        elseif l_opt_bl(i) > 0 && h_opt_bl(i) == 0 && hiredin_opt_bl(i) > 0
            case_opt_bl(i) = 3;
        elseif l_opt_bl(i) > 0 && h_opt_bl(i) == 0 && hiredin_opt_bl(i) == 0
            case_opt_bl(i) = 4;
        elseif l_opt_bl(i) == 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) > 0
            case_opt_bl(i) = 5;
        elseif l_opt_bl(i) == 0 && h_opt_bl(i) > 0 && hiredin_opt_bl(i) == 0
            case_opt_bl(i) = 6;
        elseif l_opt_bl(i) == 0 && h_opt_bl(i) == 0
            case_opt_bl(i) = 7;
        end
    end

    %----------------------------------------------
    % Estimate value of misallocation
    %----------------------------------------------

    misallocation = zeros(N,1);

    for i=1:N
        if case_opt(i) == 5 || case_opt(i) == 6   % Misallocation is zero if case 5 or 6 dominates at high SS
            misallocation(i) = 0;
        else
            misallocation(i) = profit_opt(i) - profit_opt_bl(i);
        end
    end

    misallocation(misallocation < 0) = 0;

    misallocation_untopcoded=misallocation;
    total_misallocation_untopcoded=sum(misallocation_untopcoded);

    % Top-code top 5% outliers
    misallocation_95_pcile = prctile(misallocation,95);
    misallocation(misallocation>misallocation_95_pcile) = misallocation_95_pcile;
    total_misallocation=sum(misallocation);

    %----------------------------------------------
    % Calculate cost for all HHs to reach unstable ss, i.e. overcome poverty trap (at baseline pre-transfer)
    %----------------------------------------------

    min_transfer = zeros(N,1);

    for i=1:N
        if case_opt(i) == 5 || case_opt(i) == 6   % Optimal transfer is zero if case 5 or 6 dominates at high SS
            min_transfer(i) = 0;
        else
            min_transfer(i) = uns_threshold - k_bl_pre_trans(i);
        end
    end

    min_transfer(min_transfer < 0) = 0;

    total_transfer=sum(min_transfer);

    %----------------------------------------------
    % Export results to Excel
    %----------------------------------------------

    filename1=strcat('../../Data/structural/intermediate/counterfactuals/Reducing A to point of misallocation equals cost/Simulate_reducing_A_by_pc',num2str(j));
    save(filename1);
    % tabulate(case_opt)

    %----------------------------------------------
    % Stop loop if value of misallocation = cost of eliminating trap
    %----------------------------------------------

    if total_misallocation<=total_transfer
        break
    end

end
